Глава 13 Моделирование геополей
В лекции рассмотрены детерминистические методы восстановления непрерывных поверхностей по данным в точках. Детерминистические методы интерполяции производят интерполяцию значений на основе заданной аналитической зависимости между значением в точке и значениями в пунктах с исходными данными. Это отличает их от методов геостатистических, где эта зависимость находится статистическим путем. Геостатистические методы будут рассмотрены далее.
13.1 Введение
Интерполяция в общем случае — это способ нахождения промежуточных значений величины по имеющемуся дискретному набору известных значений. В географии обычно имеют дело с двумерным случаем интерполяции — когда измерения проведены в некоторых географических локациях, и по ним нужно восстановить непрерывную картину поля распределения величины. В общем случае неизвестно, как ведет себя исследуемое явление между точками, поэтому существует бесчисленное множество вариантов интерполяции.
Методы которые производят интерполяцию на основе заданной аналитической зависимости, называют детерминистическими. Параметры этой зависимости могут быть как априори заданы пользователем, так и определяться автоматически одним из методов оптимизации — в частности, по методу наименьших квадратов.
Например, мы можем сказать, что между соседними точками показатель меняется линейным образом (здесь нужно еще указать, что мы понимаем под соседством). Такие методы достаточно просты в использовании и интерпретации. В то же время, они не учитывают статистических особенностей поведения величины между точками, которое определяется ее автокорреляционными свойствами. Методы интерполяции, которые учитывают пространственную автокорреляцию, называют геостатистическими. Они более сложны в использовании, но потенциально могут дать более достоверные результаты.
В этом модуле мы познакомимся со следующими детерминистическими методами интерполяции.
- Метод ближайшего соседа (nearest neighbor)
- Метод интерполяции на основе триангуляции
- Метод обратно взвешенных расстояний (ОВР)
- Метод радиальных базисных функций (РБФ)
- Метод иерархических базисных сплайнов (ИБС)
Интерполяцию будем рассматривать на примере данных по количеству осадков на метеостанциях в северной Италии (долина реки По и окружающие горы). Станции распределены в пространстве нерегулярно, что позволит визуально оценить чувствительность методов к этому фактору.
Прежде чем исследовать распределение показателя, необходимо проанализировать географический контекст. В этой части модуля мы воспользуемся уже известными функциями чтобы создать карту с общегеографической основой и нанести на нее пункты метеонаблюдений.
library(sf)
library(stars)
library(tmap)
library(raster)
library(plotly)
library(mapview)
library(tidyverse)
library(ggrepel)
library(dismo) # библиотека species distribution modelling
library(akima) # библиотека для интерполяции на основе триангуляции
library(gstat) # библиотека для геостатистической интерполяции, построения трендов и IDW
library(deldir) # библиотека для построения триангуляции Делоне и диаграммы Вороного
library(fields) # радиальные базисные функции (сплайны минимальной кривизны)
library(MBA) # иерархические базисные сплайны
# Убираем экспоненциальное представление больших чисел
options(scipen=999)
# Читаем слои картографической основы
cities = st_read("data/Italy_Cities.gpkg") %>% # Города
bind_cols(st_coordinates(.) %>% as_tibble())
## Reading layer `Italy_Cities' from data source `/Users/tsamsonov/GitHub/r-geo-course/data/Italy_Cities.gpkg' using driver `GPKG'
## Simple feature collection with 8 features and 37 fields
## geometry type: POINT
## dimension: XY
## bbox: xmin: 368910.4 ymin: 4930119 xmax: 686026 ymax: 5115936
## epsg (SRID): 32632
## proj4string: +proj=utm +zone=32 +datum=WGS84 +units=m +no_defs
rivers = st_read("data/Italy_Rivers.gpkg") # Реки
## Reading layer `Italy_Rivers' from data source `/Users/tsamsonov/GitHub/r-geo-course/data/Italy_Rivers.gpkg' using driver `GPKG'
## Simple feature collection with 23 features and 10 fields
## geometry type: MULTILINESTRING
## dimension: XY
## bbox: xmin: 332178.5 ymin: 4922880 xmax: 758033.3 ymax: 5121416
## epsg (SRID): 32632
## proj4string: +proj=utm +zone=32 +datum=WGS84 +units=m +no_defs
lakes = st_read("data/Italy_Lakes.gpkg") # Озера
## Reading layer `Italy_Lakes' from data source `/Users/tsamsonov/GitHub/r-geo-course/data/Italy_Lakes.gpkg' using driver `GPKG'
## Simple feature collection with 6 features and 13 fields
## geometry type: POLYGON
## dimension: XY
## bbox: xmin: 460686.2 ymin: 4938750 xmax: 757443.7 ymax: 5113557
## epsg (SRID): 32632
## proj4string: +proj=utm +zone=32 +datum=WGS84 +units=m +no_defs
# Читаем ЦМР — цифровую модель рельефа на регулярной сетке
dem = read_stars("data/gtopo.tif")
# Читаем данные об осадках
pts = read_table2("data/Rainfall.dat") %>%
st_as_sf(coords = c('x', 'y'),
crs = st_crs(cities),
remove = FALSE)
# Координаты пригодятся нам в дальнейшем
coords = st_coordinates(pts)Построение карты
# Цветовая шкала для рельефа
dem_colors = colorRampPalette(c("darkolivegreen4", "lightyellow", "orange", "firebrick"))
# Шкала высот для рельефа
dem_levels = c(0, 50, 100, 200, 500, 1000, 1500,
2000, 2500, 3000, 3500, 4000, 5000)
dem_ncolors = length(dem_levels) - 1
dem_contours = st_contour(dem, breaks = dem_levels, contour_lines = TRUE)
old = theme_set(theme_bw())
# Смотрим как выглядит результат
ggplot() +
geom_stars(data = cut(dem, breaks = dem_levels)) +
coord_sf(crs = st_crs(dem)) +
scale_fill_manual(name = 'м',
values = dem_colors(dem_ncolors),
labels = dem_levels,
na.translate = FALSE,
guide = guide_legend(label.vjust = -0.3, reverse = TRUE)) +
geom_sf(data = dem_contours, color = 'black', size = 0.2) +
geom_sf(data = rivers, color = 'midnightblue', size = 0.2) +
geom_sf(data = lakes, color = 'midnightblue', fill = 'lightblue', size = 0.2) +
geom_sf(data = pts, color = 'black', size = 0.5) +
geom_sf(data = cities, shape = 21, fill = 'white') +
geom_text_repel(data = cities, mapping = aes(x = X, y = Y, label = name), box.padding = 0.2)
13.2 Построение сетки
Любопытным свойством пакетов R, отвечающих за интерполяцию данных, является их индифферентность относительно того, в каких точках эта интерполяция будет производиться. Это может быть как регулярная растровая сетка, так и множество точек в совершенно произвольных конфигурациях. Подобная гибкость делает процесс интерполяции данных чуть более сложным, чем в ГИС-пакетах, однако способствует полному и глубокому пониманию происходящего. Вы своими руками задаете пункты, в которых следует интерполировать значения.
Для построения сетки воспользуемся функцией st_as_stars() из пакета stars, передав ей ограничивающий прямоугольник исходного множества точек и разрешение сетки:
# ПОСТРОЕНИЕ СЕТКИ ДЛЯ ИНТЕРПОЛЯЦИИ
# получим ограничивающий прямоугольник вокруг точек:
box = st_bbox(pts)
envelope = box[c(1,3,2,4)]
px_grid = st_as_stars(box, dx = 10000, dy = 10000)
ggplot() +
geom_sf(data = pts, color = 'red') +
geom_sf(data = st_as_sf(px_grid), size = 0.5, fill = NA)
Получившееся представление можно назвать сеточной моделью. В пределах каждой ячейки величина будет считаться постоянной. Ее значение будет интерполировано в центре пиксела. Можно видеть, что интерполяционная сетка слегка выходит за пределы исходного охвата множества точек. Это связано с тем, что размеры прямоугольника, ограничивающего множество точек, не кратны выбранному разрешению растра (\(10 000\) м). Проблема, однако, не столько критична, если вы выбираете достаточно подробное (малое) разрешение растра, что мы и сделаем. Зададим его равным \(1 000\) м:
13.3 Интерполяционные методы
# Цветовая шкала для осадков
rain_colors = colorRampPalette(c("white", "dodgerblue", "dodgerblue4"))
# Шкала количества осадков и соответствющее число цветов
rain_levels = seq(0, 80, by=10)
rain_ncolors = length(rain_levels)-1
rain_legend = scale_fill_manual(name = 'мм',
values = rain_colors(rain_ncolors),
guide = guide_legend(label.vjust = -0.3, reverse = TRUE, title.position = "bottom"),
labels = rain_levels,
na.value = 'white',
drop = FALSE)
rain_mapping = aes(fill = cut(rain_24, breaks = rain_levels))13.3.1 Метод ближайшего соседа (nearest neighbour)
Данный метод является простейшим по сути подходом к интерполяции. В его основе лежит построение диаграммы Вороного исходного множества точек. Считается, что в пределах каждой ячейки диаграммы значение показателя постоянно и равно значению в центре ячейки. Далее поверх диаграммы накладывается сетка интерполируемых точек и снимаются соответствующие значения:
# МЕТОД БЛИЖАЙШЕГО СОСЕДА (NEAREST NEIGHBOR)
# Диаграмма Вороного
voronoi_sf = voronoi(coords, envelope) %>%
st_as_sf() %>%
st_set_crs(st_crs(pts)) %>%
st_join(pts)
# Триангуляция Делоне
edges = pts %>%
st_union() %>%
st_triangulate()
# Визуализация
ggplot() +
geom_sf(data = voronoi_sf,
mapping = rain_mapping,
size = 0.2) +
rain_legend +
geom_sf(data = pts, color = 'red', size = 0.5) +
geom_sf(data = edges, color = 'red', size = 0.1, fill = NA)
Если есть задача конвертировать это в растр, то надо сформировать новый растр и “перенести” на него информацию с полигонов:
# Создаем растр
rnn = st_rasterize(voronoi_sf["rain_24"], px_grid)
# Визуализируем:
ggplot() +
geom_stars(data = cut(rnn, breaks = rain_levels)) +
rain_legend +
coord_sf(crs = st_crs(pts))
Видно, что полученная поверхность уже пиксельная. Для наглядности визуализируем ее в трехмерном виде. Для этого используем замечательный пакет plotly, предоставляющий интерфейс к одноименной библиотеке. Функция plot_ly(), отвечающая за построение графиков в этом пакете, требует для визуализации поверхности предоставить три компоненты:
x- вектор координат ячеек по оси \(Х\)y- вектор координат ячеек по оси \(Y\)z- матрицу значений, имеющую размеры \(length(x) \times length(y)\)
Поверхность будет раскрашиваться в различные цвета в зависимости от значений z, для управления цветами можно определить параметр colors, который должен иметь тип colorRamp:
rain_colors3d = colorRamp(c("white", "dodgerblue", "dodgerblue4"))
x = coords_grid[,'x'] %>% unique() # Получим координаты столбцов
y = coords_grid[,'y'] %>% unique() # Получим координаты строк
p = plot_ly(x = x,
y = y,
z = rnn$rain_24,
type = "surface",
colors = rain_colors3d)
layout(p, scene = list(aspectratio = list(x = 1, y = 1, z = 0.3)
))Понятное дело, что такая ступенчатая форма поверхности вряд ли соответствует реальному распределению показателя, который меняется в пространстве непрерывным образом. Можно сказать, что это первое приближение пространственного распределения.
13.3.2 Интерполяция на основе триангуляции
Интерполяция на основе триангуляции — метод интерполяции, результатом которого является уже не ступенчатая поверхность, а аппроксимированная треугольными гранями — наподобие того как объекты представлены в системах трехмерного моделирования и компьютерных играх. Триангуляция представляет собой поверхность, склеенную из треугольников, соединяющих исходные точки. Каждый треугольник является участком наклонной плоскости.
Границей триангуляции является выпуклая оболочка множества точек
Для выполнения интерполяции в произвольно заданной точке \((x, y)\) необходимо найти уравнение плоскости того треугольника который включает данную точку. В общем виде уравнение плоскости содержит четыре неизвестных коэффициента: \[ Ax + By + Cz + D = 0, \] Имея три точки \(p_1\), \(p_2\) и \(p_3\), искомые коэффициенты можно получить путем решения уравнения, левая часть которого задана в форме определителя: \[ \begin{vmatrix} x - x_1 & y - y_1 & z - z_1 \\ x_2 - x_1 & y_2 - y_1 & z_2 - z_1 \\ x_3 - x_1 & y_3 - y_1 & z_3 - z_1 \end{vmatrix} = 0 \] Коэффициенты \(A\), \(B\), \(C\) и \(D\) вычисляются заранее для каждого треугольника и хранятся вместе с триангуляцией. Получив коэффициенты нужного треугольника, искомую величину \(z(x, y)\) можно найти, выразив ее из вышеприведенного уравнения плоскости: \[ z(x, y) = -\frac{1}{C}(Ax+By+D) \]
Рассмотрим применение данного метода на практике.
# ИНТЕРПОЛЯЦИЯ НА ОСНОВЕ ТРИАНГУЛЯЦИИ (TRIANGULATION)
# Интерполируем. Параметр linear говорит о том, что показатель будет меняться вдоль ребер триангуляции линейно:
px_grid = px_grid %>%
mutate(z_linear = interpp(x = coords[,1],
y = coords[,2],
z = pts$rain_24,
xo = coords_grid[,1],
yo = coords_grid[,2],
linear = TRUE)$z)
cont_linear = st_contour(px_grid['z_linear'], breaks = rain_levels, contour_lines = TRUE)
# Смотрим как выглядит результат
ggplot() +
geom_stars(data = cut(px_grid['z_linear'], breaks = rain_levels)) +
rain_legend +
coord_sf(crs = st_crs(pts)) +
geom_sf(data = cont_linear, color = 'black', size = 0.2) +
geom_sf(data = pts, color = 'red', size = 0.5) +
geom_sf(data = edges, color = 'red', size = 0.1, fill = NA)
Обратите внимание, что все изломы (повороты) изолиний происходят на ребрах триангуляции, а внутри треугольников изолинии проходят параллельно друг другу. Каждый треугольник представляет собой фрагмент наклонной плоскости. Такой метод интерполяции, по сути, является самым простым и “честным” подходом, который близок к тому как горизонтали интерполируются вручную.
Рассмотрим поверхность в 3D:
p = plot_ly(x = x,
y = y,
z = px_grid$z_linear,
type = "surface",
colors = rain_colors3d)
layout(p, scene = list(aspectratio = list(x = 1, y = 1, z = 0.3)
))Линейная интерполяция на треугольниках, как можно видеть, выглядит достаточно угловато, хотя и существенно более правдоподобна, нежели ступенчатая поверхность, полученная методом ближайшего соседа.
Более гладкий результат можно получить, используя не линейный, а бикубический метод интерполяции на треугольниках, известный так же как метод Акимы (1978). В данном методе интерполяция базируется на следующих трёх ограничениях:
Значение функции в точке \((x, y)\) интерполируется двумерным полиномом \(5\)-й степени, который содержит \(21\) коэффициент: \[ z(x, y) = \sum_{j=0}^5 \sum_{k=0}^{5-j} q_{jk} x^j y^k \]
Значения функции, а также ее первых и вторых частных производных (\(x\), \(z_x\), \(z_y\), \(z_{xx}\), \(z_{xx}\) и \(z_{yy}\)) задаются в каждой вершине треугольника, что дает \(18\) дополнительных условий.
Частная производная функции по направлению, перпендикулярному каждой из сторон треугольника, представляет собой полином степени не более \(3\). Поскольку треугольник имеет три стороны, это ограничение порождает также еще три дополнительных условия, что в совокупности позволяет вычислить 21 искомый коэффициент.
Последнее ограничение позволяет гарантировать гладкость функции на ребрах треугольников — свойство, которое не обеспечивается обычной линейной интерполяцией.
Для того чтобы использовать методы Акимы, необходимо в функции interpp указать параметр linear = FALSE. Помимо этого, параметр extrap = TRUE говорит о том, что можно производить экстраполяцию за пределами выпуклой оболочки точек (такая возможность недоступна в линейном случае):
px_grid = px_grid %>%
mutate(z_spline = interpp(x = coords[,1],
y = coords[,2],
z = pts$rain_24,
xo = coords_grid[,1],
yo = coords_grid[,2],
linear = FALSE,
extrap = TRUE)$z)
cont_spline = st_contour(px_grid['z_spline'],
breaks = rain_levels,
contour_lines = TRUE)
# Смотрим как выглядит результат
ggplot() +
geom_stars(data = cut(px_grid['z_spline'], breaks = rain_levels)) +
rain_legend +
coord_sf(crs = st_crs(pts)) +
geom_sf(data = cont_spline, color = 'black', size = 0.2) +
geom_sf(data = pts, color = 'red', size = 0.5) +
geom_sf(data = edges, color = 'red', size = 0.1, fill = NA)
Полученные изолинии отличаются более плавным и естественным рисунком. Тем не менее, использование триангуляции все еще заметно по фестончатым изгибам изолиний на ребрах.
Смотрим наглядное представление поверхности в 3D:
p = plot_ly(x = x,
y = y,
z = px_grid$z_spline,
type = "surface",
colors = rain_colors3d)
layout(p, scene = list(aspectratio = list(x = 1, y = 1, z = 0.3)
))Можно видеть, что в данном случае получена уже гладкая поверхность, плотно натянутая на ребра триангуляции.
13.3.3 Метод обратно взвешенных расстояний (IDW)
В методе обратно взвешенных расстояний значение показателя в произвольной точке получается как средневзвешенная сумма значений в исходных точках. Веса определяются обратно пропорционально расстоянию: чем дальше исходная точка удалена, тем меньший вес она будет иметь в оценке. Формально значение функции в точке определяет согласно следующей формуле: \[ z(\mathbf{p}) = \begin{cases} \dfrac{\sum_{i = 1}^{N}{ w_i(\mathbf{p}) z_i } }{ \sum_{i = 1}^{N}{ w_i(\mathbf{p}) } }, & \text{если } d(\mathbf{p},\mathbf{p}_i) \neq 0 \text{ для всех } i, \\ z_i, & \text{если } d(\mathbf{p},\mathbf{p}_i) = 0 \text{ хотя бы для одного } i, \end{cases} \] где \(w_i(\mathbf{x}) = | \mathbf x - \mathbf x_i | ^{-p}\) — весовая функция.
Метод Шепарда — одна из наиболее распространенных модиификаций метода IDW. Веса вычисляются по формуле:
\[w_i(\mathbf{p}) = d^{-2}_i / \sum_{j=1}^n d^{-2}_j\]
Метод реализуется в R с помощью функции idw() из пакета gstat. Основным параметром метода является степень idp =, которая указывает, насколько быстро в зависимости от расстояния будет убывать вес исходной точки. По умолчанию idp = 2. При больших значениях степени (3, 4, 5, …) поверхность становится более платообразной, при меньших — островершинной.
Функция idw() принимает 4 параметра:
- Формула, указывающая название зависимой переменной и независимых переменных
- Исходные точки
- Результирующие точки
- Степень весовой функции
idp
Формулы полезны в тех случаях, когда известно (или делается предположение), что исследуемый показатель функционально связан с другой величиной. В этом случае запись Z ~ x означает, что сначала будет построена линейная регрессия \(Z(x)\) и на основе нее получена грубая оценка показателя в каждой результирующей точке. Интерполяции же будут подвергаться случайные остатки между исходными величинами в точках и теми, что получены по регрессии. Эти остатки добавляются в результирующих точках к оценке, полученной по регрессии.
С этой техникой мы познакомимся подробнее в следующем модуле при рассмотрении универсального кригинга. А пока что мы воспользуемся стандартной записью вида Z ~ 1, которая означает, что интерполироваться будет непосредственно исходная величина. В качестве Z надо указать название столбца, содержащего значения показателя. Этот столбец должен находиться в слое с исходными точками, который передается в параметр locations =. Сетка новых точек передается в параметр newdata.
Рассмотрим, как меняется вид поверхности при разных значениях idp.
# МЕТОД ОБРАТНО ВЗЕШЕННЫХ РАССТОЯНИЙ (IDW --- INVERSE DISTANCE WEIGHTED)
# Интерполируем количество осадков:
px_grid = px_grid %>%
mutate(z_idw2 = gstat::idw(rain_24 ~ 1, locations = pts, newdata = px_grid, idp = 2.0)$var1.pred,
z_idw3 = gstat::idw(rain_24 ~ 1, locations = pts, newdata = px_grid, idp = 3.0)$var1.pred,
z_idw4 = gstat::idw(rain_24 ~ 1, locations = pts, newdata = px_grid, idp = 4.0)$var1.pred,
z_idw5 = gstat::idw(rain_24 ~ 1, locations = pts, newdata = px_grid, idp = 5.0)$var1.pred)
cont_idw2 = st_contour(px_grid['z_idw2'],
breaks = rain_levels,
contour_lines = TRUE)
cont_idw3 = st_contour(px_grid['z_idw3'],
breaks = rain_levels,
contour_lines = TRUE)
cont_idw4 = st_contour(px_grid['z_idw4'],
breaks = rain_levels,
contour_lines = TRUE)
cont_idw5 = st_contour(px_grid['z_idw5'],
breaks = rain_levels,
contour_lines = TRUE)
ggplot() +
geom_stars(data = cut(px_grid['z_idw2'], breaks = rain_levels)) +
rain_legend +
coord_sf(crs = st_crs(pts)) +
geom_sf(data = cont_idw2, color = 'black', size = 0.2) +
geom_sf(data = pts, color = 'red', size = 0.5)
ggplot() +
geom_stars(data = cut(px_grid['z_idw3'], breaks = rain_levels)) +
rain_legend +
coord_sf(crs = st_crs(pts)) +
geom_sf(data = cont_idw3, color = 'black', size = 0.2) +
geom_sf(data = pts, color = 'red', size = 0.5)
ggplot() +
geom_stars(data = cut(px_grid['z_idw4'], breaks = rain_levels)) +
rain_legend +
coord_sf(crs = st_crs(pts)) +
geom_sf(data = cont_idw4, color = 'black', size = 0.2) +
geom_sf(data = pts, color = 'red', size = 0.5)
ggplot() +
geom_stars(data = cut(px_grid['z_idw5'], breaks = rain_levels)) +
rain_legend +
coord_sf(crs = st_crs(pts)) +
geom_sf(data = cont_idw5, color = 'black', size = 0.2) +
geom_sf(data = pts, color = 'red', size = 0.5)
Поскольку производная функции ОВР в каждой точке исходных данных равняется нулю (это следует из определения функции), данный метод формирует хорошо заметные замкнутые изолинии, оконтуривающие точки исходных данных. Данный эффект носит название эффекта “бычьих глаз” и относится наиболее критикуемым недостаткам метода ОВР. При увеличении степени весовой функции, однако, происходит размывание данного эффекта, выраженные в изолиниях вершины и впадины приобретают платообразный характер.
Наглядное представление о характере поверхности, получаемой методом ОВР, дает трехмерная визуализация:
p = plot_ly(x = x,
y = y,
z = px_grid$z_idw3,
type = "surface",
colors = rain_colors3d)
layout(p, scene = list(aspectratio =
list(x = 1, y = 1, z = 0.3)
))Любопытным фактом является то, что при стремлении параметра idp к плюс-бесконечности получаемая поверхность становится все более похожей не результат интерполяции методом ближайшего соседа. А именно этот метод, как мы помним, дает ступенчатую платообразную поверхность. Данный эффект легко проверить на практике, задав достаточно большой параметр idp, например \(30\):
px_grid = px_grid %>%
mutate(z_idw30 = gstat::idw(rain_24 ~ 1, locations = pts, newdata = px_grid, idp = 30.0)$var1.pred)
## [inverse distance weighted interpolation]
cont_idw30 = st_contour(px_grid['z_idw30'],
breaks = rain_levels,
contour_lines = TRUE)
ggplot() +
geom_stars(data = cut(px_grid['z_idw30'], breaks = rain_levels)) +
rain_legend +
coord_sf(crs = st_crs(pts)) +
geom_sf(data = cont_idw30, color = 'black', size = 0.2) +
geom_sf(data = pts, color = 'red', size = 0.5) +
geom_sf(data = voronoi_sf, color = 'violet', fill = NA, size = 0.2)
Рассмотрим полученную поверхность в 3D:
13.3.4 Метод радиальных базисных функций (РБФ)
Радиальной функцией называется вещественнозначная функция, значение которой зависит только от расстояния между аргументом \(\textbf{p}\) и некой фиксированной точкой в пространстве \(\textbf{c}\):
\[\phi(\textbf{p}, \textbf{c}) = \phi\big(\lVert \textbf{p} - \textbf{c}\rVert\big)\] Распространенный вариант двумерной радиальной функции возникает при \(\textbf{c} = (0, 0)\). В этом случае функция зависит только от расстояния до начала координат: \(\phi(\textbf{p}) = \phi\big(\lVert \textbf{p}\rVert\big)\)
Базисными радиальные функции называются в том случае, когда их линейная комбинация используется для аппроксимации произвольных вещественнозначных функций:
\[Z(\textbf{p}) = \sum_{i=1}^n \lambda_i \phi\big(\lVert \textbf{p} - \textbf{p}_i\rVert\big)\] Именно эта идея и используется в методе интерполяции на основе радиальных базисных функций (РБФ).
Используя условие \(Z(\textbf{p}_i) = z_i\) равенства интерполированного значения исходному, получаем систему линейных уравнений следующего вида, из которой находятся искомые коэффициенты \(\lambda_i\):
\[\begin{pmatrix} z_1 \\ \vdots\\ z_n \end{pmatrix} = \begin{pmatrix} A_{11} & \dots & A_{1n} \\ \vdots & \ddots & \vdots \\ A_{n1} & \dots & A_{nn} \end{pmatrix} \begin{pmatrix} \lambda_1 \\ \vdots\\ \lambda_n \end{pmatrix},\]
где \(A_{ij} = \phi\big(\lVert \textbf{p}_i - \textbf{p}_j\rVert\big)\), \(i, j = 1, 2, ..., n\).
Записав данную систему в матричном виде \(\mathbf{z} = \mathbf{A} \boldsymbol{\lambda}\), найти искомые коэффициенты можно путем обращения матрицы \(\mathbf{A}\) (в случае если ее определитель не равен нулю):
\(\boldsymbol{\lambda} = \mathbf{A}^{-1} \mathbf{z}\)
Метод РБФ является одним из самых гибких благодаря широким возможностям выбора радиальной функции. Недостатком же его является то, что поверхность может выходить за пределы исходного диапазона значений (хотя и обязательно проходит через исходные точки).
К числу широко используемых радиальных функций относятся:
Мультиквадрики: \(\phi(r) = \sqrt{r^2 + \delta^2}\)
Обратные мультиквадрики: \(\phi(r) = 1 / \sqrt{r^2 + \delta^2}\)
Мульти-логарифмическая: \(\phi(r) = \ln(r^2 + \delta^2)\)
Сплайны минимальной кривизны: \(\phi(r) = r^2 \ln(r^2)\)
Особого внимания среди радиальных функций заслуживают сплайны - функции, выполняющие некоторое дополнительное условие (условия) при одновременном выполнении условий интерполяции (прохождение через исходные точки). В частности, приведенный выше сплайн минимальной кривизны (thin plate spline — TPS) дает поверхность, обладающую максимально низкой кривизной между исходными точками. Это не означает что поверхность плотно натянута на исходные точки (как в триангуляции). Скорее, в ней будут отсутствовать резкие скачки и понижения, что мы видели на поверхности, построенной методом ОВР.
На языке R сплайны минимальной кривизны реализованы в пакете fields. Сначала необходимо инициализировать процесс интерполяции с помощью функции Tps(), передав ей координаты исходных точек и значения показателя в них. Дополнительно при необходимости указывается параметр scale.type = 'unscaled', который означает, что не следует масштабировать координаты исходных точек так чтобы область определения стала квадратной:
# РАДИАЛЬНЫЕ БАЗИСНЫЕ ФУНКЦИИ (RADIAL BASIS FUNCTIONS)
pred = Tps(coords, pts$rain_24, scale.type = 'unscaled')
# После этого можно интерполировать значения с помощью функции predict():
px_grid = px_grid %>%
mutate(z_tps = predict(pred, coords_grid))
# Придется расширить шкалу, так как сплайновая поверхность выходит за пределы исходных значений:
tps_breaks = seq(-10,90,by=10)
tps_ncolors = length(tps_breaks) - 1
cont_tps = st_contour(px_grid['z_tps'],
breaks = tps_breaks,
contour_lines = TRUE)
# Виузализируем результат:
ggplot() +
geom_stars(data = cut(px_grid['z_tps'], breaks = tps_breaks)) +
scale_fill_manual(name = 'мм',
values = rain_colors(tps_ncolors),
labels = paste(tps_breaks[-tps_ncolors-1], '-', tps_breaks[-1])) +
coord_sf(crs = st_crs(pts)) +
geom_sf(data = cont_tps, color = 'black', size = 0.2) +
geom_sf(data = pts, color = 'red', size = 0.5)
Можно видеть, что по плавному характеру изолиний и отсутствию артефактов в виде “бычьих глаз” интерполяция методом РБФ существенно ближе к ожидаемому распределению показателя, а также удачно сглаживает неравномерность распределения исходных данных.
Смотрим, как выглядит поверхность в 3D:
13.3.5 Метод иерархических базисных сплайнов (B-сплайнов)
В методе иерархических базисных сплайнов интерполяция прооизводится итеративно на нескольких уровнях детализации сетки: от более грубых к детальным. Задача интерполяции формулируется следующим образом.
Пусть дана область \(\Omega = \big\{(x, y)~|~0 \leq x < m,~ 0 \leq y < n \big\}\) и множество точек \(P\) нерегулярно распределенных в этой области. Для выполнения интерполяции на эту область накладывается сетка \(\Phi\) размером \((m + 3) \times (n + 3)\) — это гарантирует, что каждая точка \(P\) (в том числе, находящаяся на границе \(\Omega\)) будет окружена окрестностью из \(4 \times 4 = 16\) узлов сетки (Рис. 13.1):
Рис. 13.1: Расчетная сетка поверх области исходных данных
Пусть сетка \(\Phi\) состоит из узлов \(\phi_{ij}\) и значения индексов меняются в следующих пределах: \(i = -1, 0, ..., m + 1, ~ j = -1, 0, ..., n + 1\). Тогда функция аппроксимации в узле \(\phi_{ij}\) определяется следующим образом:
\[f(x, y) = \sum_{k=0}^3 \sum_{l=0}^3 B_k(s)B_l(t) \phi_{i+k, j+l},\] где \(i = \lfloor x \rfloor -1\), \(j = \lfloor y \rfloor -1\), \(s = x - \lfloor x \rfloor\) и \(t = y - \lfloor y \rfloor\).
Функции \(B_k\) и \(B_l\) представляют собой кубические базисные функции в форме B-сплайнов, определяемые как:
\[B_0(t) = (1-t)^3/6,\\ B_1(t) = (3t^3 - 6t^2 + 4)/6,\\ B_2(t) = (-3t^3 + 3t^2 + 3t + 1)/6,\\ B_3(t) = t^3/6,\]
где \(0 \leq t < 1\). Данные функции служат в качестве весовых коэффициентов соответствующих узлов при определении значения в точке \((x, y)\).
Интерполируемые значения \(\phi_{ij}\) находятся исходя из условия прохождения функции \(f(x, y)\) через точки исходных данных. Для значения \(z_c\) в точке \((x_c, y_c) ~|~ 1 \leq x_c, y_c < 2\) условие прохождения можно записать как:
\[z_c = \sum_{k=0}^3 \sum_{l=0}^3 w_{kl} \phi_{kl},\]
где \(w_{kl} = B_k(s) B_l(t), s = x_c - 1, t = y_c - 1\).
Поскольку существует множество значений \(\phi_{kl}\), которые удовлетворяют приведенному сооотношению, выбирается то их сочетание, которое минимизирует сумму их квадратов \(\sum_{k=0}^3 \sum_{l=0}^3 \phi_{kl}\). Стратегически это позволяет уменьшить отклонение функции \(f\) относительно нуля по области \(\Omega\).
Решая задачу по методу наименьших квадратов, можно получить следующее выражение для значения в узле сетки \(\phi_{kl}\):
\[\phi_{kl} = \frac{w_{kl} z_c}{\sum_{a=0}^3 \sum_{b=0}^3 w_{ab}}\] Функция \(f(x,y)\), построенная таким образом, будет равна \(z_c\) в точке \((x_c, y_c)\).
Приведенное выражение может быть использовано для определения значений в 16 узлах сетки, окружающих \(z_c\) в виде окрестности \(4 \times 4\). Если перейти к рассмотрению полного множества точек, окажется, что данные окрестности могут пересекаться (Рис. 13.2, слева). Поэтому при расчетах вводится понятие множества близких данных: только исходные точки, попадающие в окрестность \(4 \times 4\) относительно расчетного узла, могут оказывать влияние на его значение (Рис. 13.2, справа):
Рис. 13.2: Перекрытие окрестностей (слева) и множество близких данных (справа)
Для того, чтобы найти комппромисс между значениями, индуцируемыми в узле разными точками \(\phi_c\), минимизируется функция ошибки:
\[e(\phi_{ij}) = \sum_c (w_c \phi_{ij} - w_c \phi_c)^2.\]
Фигурирующее в ней выражение \(w_c \phi_{ij} - w_c \phi_c\) представляет собой разницу между реальным и ожидаемым вкладом \(\phi_{ij}\) в значение функции \(f\) в точке \((x_c, y_c)\). Дифференцируя \(e(\phi_{ij})\) по \(\phi_{ij}\), получаем итоговое выражение для значения в узле интерполяции \(\phi_{ij}\):
\[\phi_{ij} = \frac{\sum_c w_c^2 \phi_c}{\sum_c w_c^2}.\] Поскольку решение достигается методом наименьших квадратов, следует обратить внимание на следующий факт:
Метод иерархических базисных сплайнов является аппроксимирующим в случае если множество близких данных содержит более одной точки, и является интерполирующим в случае если множество близких данных содержит одну точку.
К достоинствам метода ИБС можно отнести:
- Поверхность получается сразу для всех узлов, нет необходимости решать систему линейных уравнений для каждого узла сетки;
- Метод является локальным: исходные точки, удаленные от текущего узла ЦМР далее чем на 2 узла, не оказывают на нее влияние. В результате этого метод ИБС получается чрезвычайно быстрым и эффективным в вычислительном плане.
- Мультимасштабность метода позволяет эффективно использовать его при интерполяции данных, распределенных кластерным образом — например, данных профилирования.
Метод иерархических базисных сплайнов доступен в пакете MBA. Чтобы использовать его, сначала необходимо подготовить исходные данные. Они должны представлять из себя матрицу из трех столбцов: X, Y, Показатель:
Метод ИБС, так же как и РБФ, предполагает по умолчанию, что область определения должна быть квадратной. Если разброс координат по осям X и Y не одинаков, поверхность будет искусственно растянута или сжата. Чтобы этого не произошло, необходимо сначала рассчитать пропорции ЦМР. Ранее мы уже создали объект envelope, который хранит крайние координаты по X и Y:
ratio = (envelope[2] - envelope[1])/(envelope[4] - envelope[3])
# После этих приготовлений можно осуществить интерполяцию
px_grid = px_grid %>%
mutate(z_bspline = mba.points(mba_data, coords_grid, n = 1, m = ratio)$xyz.est[, 'z'])
# Строим горизонтали
cont_bspline = st_contour(px_grid['z_bspline'],
breaks = tps_breaks,
contour_lines = TRUE)
# Виузализируем результат:
ggplot() +
geom_stars(data = cut(px_grid['z_bspline'], breaks = tps_breaks)) +
scale_fill_manual(name = 'мм',
values = rain_colors(tps_ncolors),
labels = paste(tps_breaks[-tps_ncolors-1], '-', tps_breaks[-1])) +
coord_sf(crs = st_crs(pts)) +
geom_sf(data = cont_bspline, color = 'black', size = 0.2) +
geom_sf(data = pts, color = 'red', size = 0.5)
Можно видеть, что метод иерархических базисных сплайнов обеспечивает некий оптимум представления поверхности. С одной стороны, он, как и метод РБФ, дает гладкую и достаточно генерализованную поверхность. С другой стороны, на участках с плотным размещением исходных данных метод ИБС раскрывает локальные нюансы поверхности, чего лишен метод РБФ, и что более типично для метода ОВР.
Наконец, рассмотрим результат в трехмерном виде:
p = plot_ly(x = x,
y = y,
z = px_grid$z_bspline,
type = "surface",
colors = rain_colors3d)
layout(p, scene = list(aspectratio =
list(x = 1, y = 1, z = 0.3)
))Итак, в настоящем модуле мы рассмотрели несколько распространенных методов детерминистической интерполяции поверхностей по данным в нерегулярно расположенных точках. В следующем модуле мы рассмотрим методы аппроксимации, которые могут быть полезны для работы с данными, обладающими высоким уровнем шума, а также для оценки пространственных трендов изменения показателя.
13.4 Аппроксимационные методы
Аппроксимационные методы используются для выявления пространственных трендов - глобальных или локальных. В зависимости от этого они и классифицируются. Полученная поверхность в каждом узле показывает средневзвешенное (типичное) значение в заданной окрестности. Таким образом, задача аппроксимации — убрать детали и выявить основные закономерности пространственного распределения. При проведении аппроксимации условие прохождения поверхности через исходные точки не применяется. В случае глобального тренда окрестность аппроксимации включает весь набор исходных точек. Локальные аппроксимации учитывают только ближайшие точки, причем общепринятым является подход, в котором окрестность определяется не расстоянием, а заданным количеством ближайших точек (или их долей от общего числа). В этом случае в области сгущения исходных данных локальная аппроксимация будет строиться по меньшей окрестности, что позволит отразить нюансы изменения показателя. И в локальных и в глобальных аппроксимациях используются обычно полиномиальные поверхности степени от \(0\) до \(3\). Коэффициенты полиномов подбираются методом наименьших квадратов для минимизации отклонения поверхности от исходных точек в заданной окрестности. В случае если степень равна \(0\), поверхность представляет из себя константу, или горизонтальную плоскость. Для степени \(1\) возможно построение наклонной плоскости. Степени \(2\) и \(3\) соответствуют квадратичным и кубическим поверхностям. Степени большего порядка для построения трендов, как правило, не используются.
13.4.1 Глобальный тренд
Построение поверхности глобального тренда можно осуществить с помощью геостатистического пакета gstat, с которым мы познакомимся в следующем модуле. Для этого необходимо сначала создать объект gstat, используя формулу (см. метод ОВР), исходные точки и степень аппроксимации. После этого аппроксимация осуществляется с помощью функции predict(). Дальнейшие действия совпадают со стандартным алгоритмом, который мы использовали ранее.
# ГЛОБАЛЬНАЯ АППРОКСИМАЦИЯ (GLOBAL APPROXIMATION)
# Создаем объект gstat и интерполируем на его основе. Столбец, указываемый в параметре formula, должен содержаться
# в наборе данных, который передается в параметр data:
px_grid = px_grid %>%
mutate(z_trend1 = predict(gstat(formula = rain_24 ~ 1, data = pts, degree = 1),
newdata = px_grid)$var1.pred,
z_trend2 = predict(gstat(formula = rain_24 ~ 1, data = pts, degree = 2),
newdata = px_grid)$var1.pred,
z_trend3 = predict(gstat(formula = rain_24 ~ 1, data = pts, degree = 3),
newdata = px_grid)$var1.pred)
# Строим горизонтали
cont_trend1 = st_contour(px_grid['z_trend1'],
breaks = rain_levels,
contour_lines = TRUE)
cont_trend2 = st_contour(px_grid['z_trend2'],
breaks = rain_levels,
contour_lines = TRUE)
cont_trend3 = st_contour(px_grid['z_trend3'],
breaks = rain_levels,
contour_lines = TRUE)
ggplot() +
geom_stars(data = cut(px_grid['z_trend1'], breaks = rain_levels)) +
rain_legend +
coord_sf(crs = st_crs(pts)) +
geom_sf(data = cont_trend1, color = 'black', size = 0.2) +
geom_sf(data = pts, color = 'red', size = 0.5)
ggplot() +
geom_stars(data = cut(px_grid['z_trend2'], breaks = rain_levels)) +
rain_legend +
coord_sf(crs = st_crs(pts)) +
geom_sf(data = cont_trend2, color = 'black', size = 0.2) +
geom_sf(data = pts, color = 'red', size = 0.5)
ggplot() +
geom_stars(data = cut(px_grid['z_trend3'], breaks = rain_levels)) +
rain_legend +
coord_sf(crs = st_crs(pts)) +
geom_sf(data = cont_trend3, color = 'black', size = 0.2) +
geom_sf(data = pts, color = 'red', size = 0.5)
Наконец, рассмотрим полученные поверхности в трехмерном виде:
p = plot_ly(x = x, y = y,
z = px_grid$z_trend1,
type = "surface",
colors = rain_colors3d)
layout(p, scene = list(aspectratio = list(x = 1, y = 1, z = 0.3)))13.4.2 Локальный тренд
Метод построения локальной регрессии изначально был разработан для построения кривых регрессии в случае когда зависимость между переменными ведет себя сложным образом и не может быть описана в терминах традиционной линейной и нелинейной регрессии — глобальных методов. В этом случае область значений независимой переменной \(X\) можно покрыть конечным числом отрезков, для каждого из которых далее находят регрессию традиционным методом — как правило, линейную или квадратичную. В классической постановке метод реализуется следующим образом. Пусть дано \(N\) точек с координатами \(X\) (независимая переменная) и \(Y\) (зависимая). Задается число \(\alpha\), которое обозначает долю от общего количества точек, которую необходимо выбрать в окрестности каждой из N точек для построения регрессии. То есть для каждой точки \(p(x)\) из исходных данных выбираются \(\alpha N\) ближайших к ней. Близость определяется как разность координат \(X\). Выбранные точки определяют окрестность \(p(x)\), в которой будет строиться локальная регрессия. Далее происходит определение параметров линейной или квадратической регрессии взвешенным методом наименьших квадратов. При использовании этого метода более близкие к \(p(x)\) точки оказывают большее влияние на коэффициенты регрессии. Построенная регрессия дает в точке \(x\) сглаженную оценку \(p'(x)\) вместо исходной \(p(x)\). Процедура повторяется для каждой из \(N\) точек. Результирующая кривая соединяет точки \(p'(x)\). При этом чем больше значение \(\alpha\), тем более сглаженный вид будет иметь кривая регрессии. Метод получил название LOWESS (Locally weighted scatterplot smoothing). В дальнейшем эта аббревиатура была редуцирована до LOESS.
В методе LOESS используются степени регрессии 0, 1, 2. Кубические и более высокие степени полиномов не применяются. При степени равной 0 метод носит название сглаживающего среднего.
Вместо координат исходных точек для построения регрессии можно использовать и произвольные координаты \(X\). В этом случае кривая будет соединять точки, полученные локальной регрессионной оценкой в заданных координатах \(X\). Именно этот принцип используется в двумерном (и многомерном) случае. Пусть даны измерения показателя в \(N\) исходных точках и задано число \(\alpha\) — сглаживающий параметр. Тогда аппроксимация показателя в каждом узле интерполяции получается путем построения поверхности тренда (см. выше) по \(\alpha N\) ближайшим исходным точкам. Как и в одномерном случае, близкие точки будут оказывать более сильное влияние на коэффициенты регрессии, чем удаленные.
Метод LOESS предоставляет широкие возможности настройки благодаря вариативности параметра сглаживания и степени регрессионного полинома.
Рассмотрим применение метода LOESS на примере данных по осадкам. Поскольку это один из базовых методов регрессионного анализа, он входит в состав базового пакета stats. Для его использования нужно вначале инициализировать параметры локальной регрессии с помощью функции loess(). Параметры задаются в следующей форме:
- Формула, содержащая названия зависимой и независимых (координаты) переменной
- Набор данных, в котором содержатся значения переменных
- Степень полинома (
degree) - Сглаживающий параметр (
span) - Необходимость нормализации координат (приведения к квадратной области определения)
# ЛОКАЛЬНАЯ АППРОКСИМАЦИЯ (LOWESS)
# 0-я степень -------------------------------------------------------------
px_grid = px_grid %>%
mutate(z_local0 = predict(loess(rain_24 ~ x + y, pts, degree = 0,
span = 0.07, normalize = FALSE), coords_grid),
z_local1 = predict(loess(rain_24 ~ x + y, pts, degree = 1,
span = 0.07, normalize = FALSE), coords_grid),
z_local2 = predict(loess(rain_24 ~ x + y, pts, degree = 2,
span = 0.07, normalize = FALSE), coords_grid))
# Визуализируем
cont_local0 = st_contour(px_grid['z_local0'],
breaks = rain_levels,
contour_lines = TRUE)
cont_local1 = st_contour(px_grid['z_local1'],
breaks = rain_levels,
contour_lines = TRUE)
cont_local2 = st_contour(px_grid['z_local2'],
breaks = rain_levels,
contour_lines = TRUE)
ggplot() +
geom_stars(data = cut(px_grid['z_local0'], breaks = rain_levels)) +
rain_legend +
coord_sf(crs = st_crs(pts)) +
geom_sf(data = cont_local0, color = 'black', size = 0.2) +
geom_sf(data = pts, color = 'red', size = 0.5)
ggplot() +
geom_stars(data = cut(px_grid['z_local1'], breaks = rain_levels)) +
rain_legend +
coord_sf(crs = st_crs(pts)) +
geom_sf(data = cont_local1, color = 'black', size = 0.2) +
geom_sf(data = pts, color = 'red', size = 0.5)
ggplot() +
geom_stars(data = cut(px_grid['z_local2'], breaks = rain_levels)) +
rain_legend +
coord_sf(crs = st_crs(pts)) +
geom_sf(data = cont_local2, color = 'black', size = 0.2) +
geom_sf(data = pts, color = 'red', size = 0.5)
Рассмотрим результаты в 3D:
p = plot_ly(x = x, y = y,
z = px_grid$z_local0,
type = "surface",
colors = rain_colors3d)
layout(p, scene = list(aspectratio = list(x = 1, y = 1, z = 0.3)))
p = plot_ly(x = x, y = y,
z = px_grid$z_local1,
type = "surface",
colors = rain_colors3d)
layout(p, scene = list(aspectratio = list(x = 1, y = 1, z = 0.3)))
p = plot_ly(x = x, y = y,
z = px_grid$z_local2,
type = "surface",
colors = rain_colors3d)
layout(p, scene = list(aspectratio = list(x = 1, y = 1, z = 0.3)))Можно заметить, что с увеличением степени полинома поверхность все более точно аппроксимирует исходные данные — там, где достаточно большое количество исходных точек. В то же время, появляются нежелательные экстраполяции в приграничных областях, слабо обеспеченных измерениями, что можно наблюдать на последнем рисунке (степень = 2) в северо-западной части. Поэтому нельзя однозначно сказать, что более высокая степень обеспечивает лучшие результаты аппроксимации. Точность аппроксимации правильно регулировать не степенью полинома, а увеличением и уменьшением сглаживающего параметра альфа.
Рассмотрим это на примере линейной аппроксимации при \(\alpha = 0.05, 0.1, 0.2\):
px_grid = px_grid %>%
mutate(z_local1_005 = predict(loess(rain_24 ~ x + y, pts, degree = 1,
span = 0.05, normalize = FALSE), coords_grid),
z_local1_01 = predict(loess(rain_24 ~ x + y, pts, degree = 1,
span = 0.1, normalize = FALSE), coords_grid),
z_local1_02 = predict(loess(rain_24 ~ x + y, pts, degree = 1,
span = 0.2, normalize = FALSE), coords_grid))
# Визуализируем
cont_local1_005 = st_contour(px_grid['z_local1_005'],
breaks = rain_levels,
contour_lines = TRUE)
cont_local1_01 = st_contour(px_grid['z_local1_01'],
breaks = rain_levels,
contour_lines = TRUE)
cont_local1_02 = st_contour(px_grid['z_local1_02'],
breaks = rain_levels,
contour_lines = TRUE)
ggplot() +
geom_stars(data = cut(px_grid['z_local1_005'], breaks = rain_levels)) +
rain_legend +
coord_sf(crs = st_crs(pts)) +
geom_sf(data = cont_local1_005, color = 'black', size = 0.2) +
geom_sf(data = pts, color = 'red', size = 0.5)
ggplot() +
geom_stars(data = cut(px_grid['z_local1_01'], breaks = rain_levels)) +
rain_legend +
coord_sf(crs = st_crs(pts)) +
geom_sf(data = cont_local1_01, color = 'black', size = 0.2) +
geom_sf(data = pts, color = 'red', size = 0.5)
ggplot() +
geom_stars(data = cut(px_grid['z_local1_02'], breaks = rain_levels)) +
rain_legend +
coord_sf(crs = st_crs(pts)) +
geom_sf(data = cont_local1_02, color = 'black', size = 0.2) +
geom_sf(data = pts, color = 'red', size = 0.5)
Сравниваем результаты в трехмерном виде:
p = plot_ly(x = x, y = y,
z = px_grid$z_local1_005,
type = "surface",
colors = rain_colors3d)
layout(p, scene = list(aspectratio = list(x = 1, y = 1, z = 0.3)))13.5 Интерполяция по ареалам
В некоторых случаях необходимо осуществить так называемую интерполяцию по ареалам. Данный метод применяется в тех случаях, когда исходная информация привязана не к точечным, а к площадным объектам. Задача заключается в том, чтобы с одной площадной сетки перенести на другую (как правило, регулярную, обладающую большей дискретностью). Необходимость подобного преобразования может быть обусловлена следующими (но и не только) причинами:
- метод анализа (например, моделирование диффузии) предполагает, что данные распределены по регулярной сетке, в то время как исходная сетка нерегулярна.
- необходимо обеспечить сравнимость пространственных распределений показателя для разных территорий, в то время как дробность исходного территориального деления существенно меняется в пространстве.
Метод интерполяции по ареалам реализуется средствами функции st_interpolate_aw() из пакета sf. Данной функции необходимо подать исходную и целевую полигональную сетку, а также указать тип параметра: интенсивный или экстенсивный:
— экстенсивные параметры суммируются и делятся при агрегировании/агрегировании территориальных единиц. Например, площадь, покрытая лесом или численность населения — это экстенсивный параметр. - интенсивные параметры осредняются или остаются постоянными при агрегировании/дизагрегировании территориальных единиц. Например, густота древостоя и плотность населения — интенсивные параметры.
Рассмотрим это метод интерполяции на примере данных по графствам Северной Каролины (показатель — количество новорожденных в 1974 году). Для расчета векторной регулярной сетки используем функцию st_make_grid() из пакета sf.
# Данные по Северной Каролине
nc = st_read(system.file("shape/nc.shp", package="sf"))
## Reading layer `nc' from data source `/Library/Frameworks/R.framework/Versions/3.6/Resources/library/sf/shape/nc.shp' using driver `ESRI Shapefile'
## Simple feature collection with 100 features and 14 fields
## geometry type: MULTIPOLYGON
## dimension: XY
## bbox: xmin: -84.32385 ymin: 33.88199 xmax: -75.45698 ymax: 36.58965
## epsg (SRID): 4267
## proj4string: +proj=longlat +datum=NAD27 +no_defs
cells = st_make_grid(nc, cellsize = 0.25)
birth = st_interpolate_aw(nc["BIR74"],
cells,
extensive = TRUE)
# исходное распределение
tm_shape(nc) +
tm_polygons('BIR74',
style = 'jenks',
palette = 'viridis') +
tm_shape(cells) +
tm_borders(col = 'white')
# пересчет на регулярную сетку
tm_shape(birth) +
tm_polygons('BIR74',
style = 'jenks',
palette = 'viridis') +
tm_shape(nc) +
tm_borders(col = 'white')
13.6 Краткий обзор
Для просмотра презентации щелкните на ней один раз левой кнопкой мыши и листайте, используя кнопки на клавиатуре:
Презентацию можно открыть в отдельном окне или вкладке браузере. Для этого щелкните по ней правой кнопкой мыши и выберите соответствующую команду.
13.7 Контрольные вопросы и упражнения
13.7.1 Вопросы
- Сформулируйте отличие интерполяционных и аппроксимационных методов восстановления поверхностей по данным в нерегулярно расположенных точках.
- Какой метод интерполяции подходит для работы с категориальными (качественными) данными?
- При каких условиях в методе интерполяции на основе триангуляции возможна оценка величины за пределами выпуклой оболочки исходного множества точек?
- Назовите основной недостаток метода обратно взвешенных расстояний.
- Каким требованиям должна удовлетворять радиальная базисная функция?
- Если дана функция, являющаяся радиальной базисной, можно ли в качестве таковой использовать обратную к ней? Как это реализуется математически?
- В каком методе интерполяции используются сплайны с натяжением?
- Сформулируйте основные принципы восстановления поверхности методом иерархических базисных сплайнов.
- Полиномы каких степеней обычно используются в задачах аппроксимации на практике?
- Каким образом реализуется метод локальной регрессии в приложении к пространственным данным?
- В чем заключается отличие экстенсивных и интенсивных пространственных переменных?
- Расскажите алгоритм, лежащий в основе пикнофилактического метода интерполяции по ареалам.
- Перечислите пакеты программной среды R, которые можно использовать для восстановления поверхностей по данным в нерегулярно расположенных точках.
- Какие функции R позволяют построить регулярную сетку? Как преобразовать эту сетку в растр?
- Изложите последовательность действий, необходимую для трехмерной визуализации поверхности средствами библиотеки plotly.
13.7.2 Упражнения
Загрузите данные дрейфующих буев ARGO на акваторию Северной Атлантики за 30 января 2010 года. Постройте поля распределения солености и температуры методами радиальных базисных функций и обратно взвешенных расстояний с размером ячейки 50 км. Визуализируйте их средствами tmap и сравните результаты. Какой метод интерполяции дает, на ваш взгляд, более правдоподобный результат?
Подсказка: перед построением сетки интерполяции странсформируйте данные в проекцию Меркатора. Чтобы полученное поле распределения покрывало только акваторию, маскируйте полученный растр с использованием слоя ocean из набора данных Natural Earth. Перед выполнением маскирования преобразуйте мультиполигон в обычные полигоны (в противном случае маскирование отработает некорректно).
Используя пакет eurostat, загрузите с портала Евростата таблицу по объему производства коровьего молока, а также соответствующие ей единицы деления 2-го уровня (NUTS 2). Используя метод интерполяции по ареалам, пересчитайте эти данные на регулярную сетку с шагом 100 км и визуализируйте средствами tmap. Сравните результаты с оригинальным распределением.
Подсказка: перед построением сетки интерполяции обрежьте данные на территорию континентальной Европы по координатам углов, используя функцию
st_crop(). Далее трансформируйте данные в коническую равнопромежуточную проекцию с параметрами+proj=eqdc +lon_0=20 +lat_1=43 +lat_2=62. После этого присоедините статистику по полю geo (метки c помощьюlabel_eurostat()можно не назначать) и произведите интерполяцию по ареалам.
| Самсонов Т.Е. Визуализация и анализ географических данных на языке R. М.: Географический факультет МГУ, 2020. DOI: 10.5281/zenodo.901911 |